#include "G4SystemOfUnits.hh"
#include "MyCsIBirksLaw.hh"
#include "PGGeneratorList.hh"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <filesystem>
#include <string>

MyCsIBirksLaw::MyCsIBirksLaw(G4int verbose_)
{
    m_verboseLevel = verbose_;
    m_E0 = 12.97;   // keV
    m_kB = 3.311*std::pow(10,-6);   //g * (keV * cm^2)^(-1)
    m_rho = 4.53;   // g/cm^3
    m_IfNeutrino = false;
    if(PGGeneratorList::GetInstance()->GetPGPSConfig().ParticleGunType == G4String("Nucleus"))
    {
        m_IfNeutrino = true;
        std::string dEdxPath = std::string(NUHUNTER_SOURCE_DIR) + std::string("/source/1_Simulation/source/PhysicsList/include/CsIBirksdEdx");
        for (const auto & entry : std::filesystem::directory_iterator(dEdxPath))
        {
            std::string path = entry.path();
            if(path.find(".txt") == std::string::npos)
            {
                continue;
            }
            std::string ParticleName = path;
            ParticleName.erase(path.find(".txt"));
            ParticleName.erase(0, path.rfind("/") + 1);

            std::vector<Double_t> Energy;   Energy.clear();
            std::vector<Double_t> dEdxVec;     dEdxVec.clear();
            std::ifstream InputCSV(path.c_str());
            std::string Line;
            while (std::getline(InputCSV, Line))
            {
                std::string Number;
                double dedx = 0.;
                std::istringstream ReadLine(Line);
                std::pair<Double_t, Double_t> EnergyVSLumi;
                std::getline(ReadLine, Number, ',');
                Energy.push_back(std::atof(Number.c_str()));
                std::getline(ReadLine, Number, ',');
                dedx = std::atof(Number.c_str());
                std::getline(ReadLine, Number, ',');
                dedx += std::atof(Number.c_str());
                dEdxVec.push_back(dedx);
            }
            m_dEdxMap[ParticleName] = new TGraph(Energy.size(), &(Energy[0]), &(dEdxVec[0]));
        }
    }

    #ifdef _MyCsIBirksLaw_DBG_
        m_outputFile.open("CsIBirksLawDBG.csv", std::ios::out);
        m_outputFile << std::resetiosflags(std::ios::scientific); 
        m_outputFile.setf(std::ios::scientific, std::ios::floatfield); m_outputFile.precision(8);
        m_outputFile << "Erec,dEdx,QF" << std::endl;
        //m_outputFile << "Erec,Edep,StepLength,QF" << std::endl;
    #endif
}

MyCsIBirksLaw::~MyCsIBirksLaw()
{
    #ifdef _MyCsIBirksLaw_DBG_
        m_outputFile.close();
    #endif
    if(m_IfNeutrino)
    {
        for(const auto& dEdx : m_dEdxMap)
        {
            delete dEdx.second;   
        }
    }
}

double MyCsIBirksLaw::QF(const G4Step& aStep)
{
    if(m_IfNeutrino == false)
    {
        //std::cout << "Warning!!! ---->MyCsIBirksLaw::QF(): not use Nucleus Particle Gun, return 1" << std::endl;
        return 1.;
    }
    double Erec = aStep.GetPreStepPoint()->GetKineticEnergy() / keV;
    //double Edep = aStep.GetTotalEnergyDeposit() / keV;
    //double StepLength = aStep.GetStepLength() / cm;
    double dedx = 0.;
    if(m_IfNeutrino == true)
    {
        const G4Track *aTrack = aStep.GetTrack();
        G4ParticleDefinition *particle = aTrack->GetDefinition();
        std::string ParticleName = particle->GetParticleName();
        if (m_dEdxMap.find(ParticleName) == m_dEdxMap.end())
        {
            std::cout << "Warning!!! ---->MyCsIBirksLaw::QF(): Particle " << ParticleName << " Not found" << std::endl;
            return 1;
        }
        else
        {
            dedx = m_dEdxMap[ParticleName]->Eval(Erec);
        }
    }

    double QF = (1 - std::exp(-Erec / m_E0)) * (1. / (m_kB * dedx));
    //double QF = (1 - std::exp(-Erec / m_E0)) * (1. / (m_kB * (Edep / StepLength / m_rho)));

    if(m_verboseLevel)
    {
        std::cout << "---->MyCsIBirksLaw::QF()" << std::endl;
        std::cout << "Erec: " << Erec << "\tdEdx: " << dedx << "\tQF:" << QF << std::endl;
        //std::cout << "Erec: " << Erec << " keV\tEdep: " << Edep << " keV\tStepLength: " << StepLength << " cm\tQF: " << QF << std::endl;
    }

    #ifdef _MyCsIBirksLaw_DBG_
        m_outputFile << Erec << "," << dedx << "," << QF << "\n";
    #endif

    return QF;
}